Waldmann (2009b) reported that his implementation is faster
and results in Markov chains with better mixing (smaller autocorrelations)
in comparison to Damgaard (2007b). The code proposed by
Waldmann (2009b) and here AnimalModelBUGSVar) first
computes the parent average and precision (1 / variance) to sample
from the normal distribution with these two parameters. Damgaard (2007b),
on the other hand, first sampled from the standard normal distribution,
scaled by the square root of a within-family segregation variance
coefficient added the parent average, and at the end scaled by the
additive genetic standard deviation. The approach of Damgaard (2007b)
is essentially the transformed (a priori independent) sampler
for additive genetic values (e.g. Waagepetersen, N., and Sorensen (2008), and Shariati and Sorensen (2009b)).
To highlight the differences between the proposed implementations,
the code by Damgaard (2007b) was rewritten in two variants
(functions Variantes 1 and Variantes 2) and compared for speed and mixing
with the code presented here AnimalModelBUGSVar). The original
implementations of Damgaard (2007b) and Waldmann (2009b)
were not tested because they are not generic.
See bellow two variants of the animal model with transformed
(a priori independent) additive genetic values
Additionally, two other parameterizations of the animal model are
interesting which involve the incorporation of NRM into the design
matrix \(\mathbf{Z}\) to remove a priori dependence between
additive genetic values, e.g., \(\mathbf{u}|\sigma_{a}^{2}\sim N\left(\mathbf{0},\mathbf{I}\sigma_{a}^{2}\right)\).
One approach is to use left Cholesky factor \(\mathbf{Z}^{*}=\mathbf{Z}\mathbf{L}=\mathbf{Z}\mathbf{T}\mathbf{W}^{\frac{1}{2}}\)
(Quaas, Anderson, and Gilmour 1984b), where original \(\mathbf{a}\) can be computed
using \(\mathbf{a}=\mathbf{T}\mathbf{W}^{\frac{1}{2}}\mathbf{u}\)
or more efficiently with solving \(\mathbf{T}^{-1}\mathbf{a}=\mathbf{W}^{\frac{1}{2}}\mathbf{u}\)
due to the sparseness of \(\mathbf{T}^{-1}\) (e.g. Waagepetersen, N., and Sorensen (2008)).
In a similar manner, Waldmann et al. (2008a) used a singular value
decomposition of NRM, i.e., \(\mathbf{A}=\mathbf{U}\mathbf{S}\mathbf{U}^{T}\),
where \(\mathbf{U}\) is a dense orthonormal matrix \(\left(\mathbf{U}\mathbf{U}^{T}=\mathbf{I}\right)\)
and \(\mathbf{S}\) is a diagonal matrix. Then the parameterization
is \(\mathbf{Z}^{*}=\mathbf{Z}\mathbf{U}\mathbf{S}^{\frac{1}{2}}\)
and \(\mathbf{a}=\mathbf{U}\mathbf{S}^{\frac{1}{2}}\mathbf{u}\).
JAGS implementation for such a parameterization AnimalModelBUGSZStarVar)
is straightforward because the additive genetic values become independent and
fitted with a homogenous variance (precision), while the pedigree
information is incorporated into the design matrix \(\mathbf{Z}^{*}\).
However, this means that the whole dense matrix \(\mathbf{Z}^{*}\)of
order \(nY\times nI\) must be passed to JAGS. Initial tests showed
that the CPU time is so long that using these two JAGS approaches is applicable only for pedigrees with a handful of members.
A similar conclusion also holds for the augmented parameterization of
the mixed model by Harville (1986a). Therefore, these approaches
were not tested for speed and mixing.
Consider now the code below, which shows an animal model with a reparameterized design matrix \(\mathbf{Z}\).
Implementations (functions AnimalModelBUGSVar, Variantes 1, and Variantes 2)
were tested using our simple example described in section 1. Here, the following
two approaches were used: one with variances assumed known and one with
variances unknown in addition to other model parameters (\(\mathbf{b}\)
and \(\mathbf{a}\)).
Sampling time reported by JAGS was collected, and for each parameter, an effective sample size (\(n_{eff}\)) was computed
using the CODA program Plummer et al. (2009). The last parameter describes
the amount of information obtained in the collected samples from MCMC
sampling. The effective sample size is almost always smaller than the
sample size due to autocorrelations between the Markov chain’s subsequent samples (e.g. Kass et al. (1998)). Here we are interested in MCMC algorithms’ efficiency for generating accurate posterior samples in a Bayesian analysis. We are using the MCMC pace or computing cost per sample, which measures the average time needed to generate each effectively independent sample, to summarise algorithm efficiency. Thus, the computing cost per \(N_{s}\) samples for the parameter \(p\) \(\left(c_{p}\right)\) was evaluated as \(c_p=N_{s} \times t_{n}/n_{eff}^{(p)}\), where \(t_{n}\) is the time spent for drawing \(n\) samples. Considering the model has \(P\) parameters, we can summarise the average computing cost per sample as \(\bar{c}=\frac{N_{s} \times t_{n}}{P}\sum_{p=1}^{P}1/n_{eff}^{(p)}\), where \(n_{eff}^{(p)}\) is the \(n_{eff}\) for the \(p\)th parameter. However, often one or few parameters may have a worse mixing than the others and, consequently, limit deciding which model should be adequately chosen. For this reason, we are also considering to include maximum computing cost per \(N_{s}\) samples (\(c_{max} = \max\lbrace{c_{p}\rbrace}_{p=1}^{P}\)) that represents the longest time taken to generate each effectively independent sample.
All computations were performed on a Dell Inspiron 17 7000 with 10\(^{\tiny \mbox{th}}\) Generation Intel\(\circledR\) Core\(^{\tiny\mbox{TM}}\) i7 processor, 1.80GHz \(\times\) four-processor speed, 16GB random access memory (RAM) plus 20GB of swap space, 64-bit integers, and the platform used is a Linux Mint 19.2 Cinnamon system version 5.2.2-050202-generic.
Example 1
We are running three chains, each of one has chain length, burn-in, and thinning rate were arbitrarily set to 20,000, 30,000, and 15, respectively, and the data example, which can be seen below:
Tables 10.1 and 10.2 show the implementation using transformed additive genetic values for known and unknown variances, respectively (Damgaard (2007b);
functions Variantes 1, and Variantes 2). In general, models m1, m2, m7, and m8 show the lowest amount of computing time while m3, m5 and m7 the lowest computing cost per sample independently if variance components are known or unknown.
Table 10.1: Effect of animal model (known variances) implementation in JAGS using the mean and standard error of sampling time (t) in seconds, computing cost per 10,000 samples (\(c_{p}=10,000 \times t_{n}/n_{eff}\)), mean computing cost per 10,000 samples (\(\bar{c}\)), and maximum computing cost per 10,000 samples (\(c_{max}\)) based on 50 models fitted using three chains each considering a length of 20,000, burn-in set to 30,000, and thinning to 15
|
|
Time
|
Computing cost per sample (c)
|
|
|
Model
|
Type
|
Mean
|
SE
|
a[1]
|
a[2]
|
a[3]
|
a[4]
|
a[5]
|
a[6]
|
a[7]
|
a[8]
|
a[9]
|
b[1]
|
b[2]
|
cost_Mean
|
cost_max
|
|
m1
|
No Inbreeding
|
1.2168
|
0.0078
|
0.2033
|
0.2281
|
0.2325
|
0.2237
|
0.2265
|
0.2346
|
0.2349
|
0.2058
|
0.2253
|
0.2175
|
0.2239
|
0.2233
|
0.2349
|
|
m2
|
Inbreeding
|
1.1781
|
0.0285
|
0.1986
|
0.2238
|
0.2300
|
0.2186
|
0.2235
|
0.2313
|
0.2342
|
0.1992
|
0.2212
|
0.2123
|
0.2189
|
0.2192
|
0.2342
|
|
m3
|
Inbreeding + Cholesky Decomposition
|
1.5529
|
0.0349
|
0.2588
|
0.2570
|
0.2568
|
0.2588
|
0.2588
|
0.2588
|
0.2631
|
0.2531
|
0.2588
|
0.2588
|
0.2599
|
0.2584
|
0.2631
|
|
m4
|
Inbreeding + Cholesky Decomposition (Variant 1)
|
2.7878
|
0.0725
|
0.4646
|
0.4614
|
0.4610
|
0.4646
|
0.4646
|
0.4646
|
0.4724
|
0.4543
|
0.4646
|
0.4646
|
0.4666
|
0.4640
|
0.4724
|
|
m5
|
Inbreeding + Cholesky Decomposition (Variant 2)
|
2.7865
|
0.0915
|
0.4644
|
0.4644
|
0.4644
|
0.4644
|
0.4644
|
0.4644
|
0.4594
|
0.4505
|
0.4703
|
0.4655
|
0.4644
|
0.4633
|
0.4703
|
|
m6
|
Inbreeding + Cholesky Decomposition + Reording Pedigree
|
1.4767
|
0.0267
|
0.2461
|
0.2481
|
0.2481
|
0.2461
|
0.2490
|
0.2416
|
0.2501
|
0.2407
|
0.2543
|
0.2461
|
0.2461
|
0.2469
|
0.2543
|
|
m7
|
Inbreeding + Cholesky Decomposition + Zstar
|
1.2577
|
0.0254
|
NA
|
0.2096
|
0.2096
|
0.2070
|
0.2096
|
0.2132
|
NA
|
NA
|
0.2096
|
0.2096
|
0.2096
|
0.2097
|
0.2132
|
|
m8
|
Inbreeding + Zstar + SVD
|
1.2598
|
0.0318
|
NA
|
0.2100
|
0.2009
|
0.2096
|
0.2100
|
0.2127
|
NA
|
NA
|
0.2100
|
0.2100
|
0.2100
|
0.2091
|
0.2127
|
Table 10.2: Effect of animal model (unknown variances) implementation in JAGS using the mean and standard error of sampling time (t) in seconds, computing cost per 10,000 samples (\(c_{p}=10,000 \times t_{n}/n_{eff}\)), mean computing cost per 10,000 samples (\(\bar{c}\)), and maximum computing cost per 10,000 samples (\(c_{max}\)) based on 50 models fitted using three chains each considering a length of 20,000, burn-in set to 30,000, and thinning to 15
|
|
Time
|
Computing cost per sample (c)
|
|
|
Model
|
Type
|
Mean
|
SE
|
a[1]
|
a[2]
|
a[3]
|
a[4]
|
a[5]
|
a[6]
|
a[7]
|
a[8]
|
a[9]
|
b[1]
|
b[2]
|
cost_Mean
|
cost_max
|
|
m1
|
No Inbreeding
|
2.2973
|
0.0847
|
1.0568
|
4.9257
|
5.0535
|
5.6073
|
5.8800
|
5.5331
|
3.7477
|
0.7853
|
5.0837
|
3.7550
|
4.3859
|
4.1649
|
5.8800
|
|
m2
|
Inbreeding
|
2.5375
|
0.0819
|
1.1750
|
5.5892
|
5.7894
|
6.3820
|
6.6601
|
6.1996
|
4.5385
|
0.9712
|
5.8833
|
4.2892
|
5.0148
|
4.7720
|
6.6601
|
|
m3
|
Inbreeding + Cholesky Decomposition
|
2.6787
|
0.0947
|
3.2560
|
4.6798
|
4.9114
|
4.7885
|
5.0266
|
4.8739
|
4.0257
|
2.3665
|
4.8109
|
3.6795
|
3.9214
|
4.2127
|
5.0266
|
|
m4
|
Inbreeding + Cholesky Decomposition (Variant 1)
|
2.8721
|
0.0823
|
3.4911
|
5.0177
|
5.2661
|
5.1343
|
5.3896
|
5.2259
|
4.3164
|
2.5374
|
5.1583
|
3.9452
|
4.2046
|
4.5170
|
5.3896
|
|
m5
|
Inbreeding + Cholesky Decomposition (Variant 2)
|
8.0390
|
0.2557
|
10.5776
|
12.9997
|
16.1426
|
13.6625
|
15.3945
|
16.1816
|
11.7857
|
9.1990
|
14.3810
|
7.2712
|
7.9202
|
12.3196
|
16.1816
|
|
m6
|
Inbreeding + Cholesky Decomposition + Reording Pedigree
|
2.3657
|
0.1066
|
3.0431
|
4.6697
|
4.6734
|
4.6468
|
4.7125
|
4.6799
|
4.0322
|
2.5079
|
4.6762
|
3.5489
|
3.6362
|
4.0752
|
4.7125
|
|
m7
|
Inbreeding + Cholesky Decomposition + Zstar
|
2.6427
|
0.0696
|
NA
|
5.4918
|
4.5485
|
4.3096
|
4.1305
|
3.6857
|
NA
|
NA
|
4.3061
|
4.2541
|
4.1676
|
4.3617
|
5.4918
|
|
m8
|
Inbreeding + Zstar + SVD
|
2.2873
|
0.0841
|
NA
|
4.8852
|
3.3503
|
3.1976
|
1.0630
|
2.9674
|
NA
|
NA
|
2.7419
|
3.3785
|
3.5778
|
3.1452
|
4.8852
|
The original implementation of Damgaard ((Damgaard 2007b); variant 2 in function Variantes 2) is the most time-consuming approach and can be even more time consuming when variances are inferred.
Although slower, the implementation with the transformed additive
genetic values constantly produced chains with a considerably larger
effective sample size due to `a prior independence among
additive genetic values. Different implementations were compared by
the computing cost per effective sample size (Tables 10.1 and 10.2)
as suggested by Waagepetersen, N., and Sorensen (2008). When variances were assumed to be known, the cost did not differ a lot between implementations, though, on average, a standard implementation was the best one. This implementation was also the best when variances were unknown.
Shariati and Sorensen (2009b) performed
extensive tests of various MCMC samplers for animal model and concluded
that none of the samplers was the best in all scenarios, but transformed
and blocked samplers could be beneficial when there are large posterior
correlations between the parameters.
Computing Cost
We are also running three chains with chain length, burn-in, and thinning rate arbitrarily set per chain to 5,000, 1,000, and 10, respectively, which give us 3,000 simulation draws in total. The results for some models are presented below, and further discussion will be addressed at the end of the section.
- Basic Animal model
# Variance unknown
AnimalModelBUGSVar
## function()
## {
## ## Priors - precision parameters
## tau2e ~ dgamma(0.001, 0.001)
## tau2a ~ dgamma(0.001, 0.001)
## sigma2e <- pow(tau2e, -1)
## sigma2a <- pow(tau2a, -1)
##
## ## Additive genetic values
## for(k in 1:nI) {
## a[id[k]] ~ dnorm(pa[id[k]], Xtau2a[id[k]])
## pa[id[k]] <- 0.5 * (a[fid[k]] + a[mid[k]])
## Xtau2a[id[k]] <- winv[id[k]] * tau2a
## }
## a[nU] <- 0 # NULL (zero) holder
##
## ## Fixed Effects
## for(j in 1:nB) { b[j] ~ dnorm(0, 1.0E-6) }
##
## ## Phenotypes
## for(i in 1:nY) {
## y[i] ~ dnorm(mu[i], tau2e)
## mu[i] <- inprod(x[i, ], b[]) + a[idy[i]]
## }
## }
- Without accounting for inbreeding:
dataNoInb <- do.call('rbind',
lapply(c("./Results/Large_Example/ESSNoInb.RDS",
"./Results/Large_Example/ESSNoInbVar.RDS"), readRDS))
dataNoInb$model <- as.factor(dataNoInb$model)
g1 <- dataNoInb %>%
ggplot(aes(y= cost_mean, x= Generation)) +
geom_point(aes(colour = par)) +
geom_line(aes(linetype = par, colour = par)) +
xlab("Generation") +
ylab("Mean computing cost per sample (seconds)")+
ylim(0, max(dataNoInb$cost_max)) +
facet_wrap(~model)+
theme_bw() +
theme(legend.position = "none")
g2 <- dataNoInb %>%
ggplot(aes(y= cost_max, x= Generation)) +
geom_point(aes(colour = par)) +
geom_line(aes(linetype = par, colour = par)) +
xlab("Generation") +
ylim(0, max(dataNoInb$cost_max)) +
ylab("Maximum computing cost per sample (seconds)")+
facet_wrap(~model)+
theme_bw()
g1+g2
- Accounting for inbreeding:
dataInb <- do.call('rbind',
lapply(c("./Results/Large_Example/ESSInb.RDS",
"./Results/Large_Example/ESSInbVar.RDS"), readRDS))
dataInb$model <- as.factor(dataInb$model)
g1 <- dataInb %>%
ggplot(aes(y= cost_mean, x= Generation)) +
geom_point(aes(colour = par)) +
geom_line(aes(linetype = par, colour = par)) +
xlab("Generation") +
ylab("Mean computing cost per sample (seconds)")+
ylim(0, max(dataInb$cost_max)) +
facet_wrap(~model)+
theme_bw() +
theme(legend.position = "none")
g2 <- dataInb %>%
ggplot(aes(y= cost_max, x= Generation)) +
geom_point(aes(colour = par)) +
geom_line(aes(linetype = par, colour = par)) +
xlab("Generation") +
ylim(0, max(dataInb$cost_max)) +
ylab("Maximum computing cost per sample (seconds)")+
facet_wrap(~model)+
theme_bw()
g1+g2
- Animal model using Cholesky decompostion of the NRM
# VAriance unknown
AnimalModelBUGSCholVar
## function()
## {
## ## Variance priors
## tau2e ~ dgamma(0.001, 0.001)
## tau2a ~ dgamma(0.001, 0.001)
## sigma2e <- 1 / tau2e
## sigma2a <- 1 / tau2a
##
## ## Additive genetic values
## for(k in 1:nI) {
## u[k] ~ dnorm(0, tau2a)
## a[id[k]] <- u[k] * wsqr[id[k]] + 0.5 * (a[fid[k]] + a[mid[k]])
## }
## a[nU] <- 0 # NULL (zero) holder
##
## ## Location priors
## for(j in 1:nB) { b[j] ~ dnorm(0, 1.0E-6) }
##
## ## Phenotypes
## for(i in 1:nY) {
## y[i] ~ dnorm(mu[i], tau2e)
## mu[i] <- inprod(x[i, ], b[]) + a[idy[i]]
## }
## }
dataChol <- do.call('rbind',
lapply(c("./Results/Large_Example/ESSChol.RDS",
"./Results/Large_Example/ESSCholVar.RDS"), readRDS))
dataChol$model <- as.factor(dataChol$model)
g1c <- dataChol %>%
ggplot(aes(y= cost_mean, x= Generation)) +
geom_point(aes(colour = par)) +
geom_line(aes(linetype = par, colour = par)) +
xlab("Generation") +
ylab("Mean computing cost per sample (seconds)")+
ylim(0, max(dataChol$cost_max)) +
facet_wrap(~model)+
theme_bw() +
theme(legend.position = "none")
g2c <- dataChol %>%
ggplot(aes(y= cost_max, x= Generation)) +
geom_point(aes(colour = par)) +
geom_line(aes(linetype = par, colour = par)) +
xlab("Generation") +
ylim(0, max(dataChol$cost_max)) +
ylab("Maximum computing cost per sample (seconds)")+
facet_wrap(~model)+
theme_bw()
g1c+g2c
- Animal model variant 1 using Cholesky decompostion of the NRM
# Variance unknown
AnimalModelBUGSVariant1Var
## function()
## {
## ## Priors - precision parameters
## tau2e ~ dgamma(0.001, 0.001)
## tau2a ~ dgamma(0.001, 0.001)
## sigma2e <- pow(tau2e, -1)
## sigma2a <- pow(tau2a, -1)
##
## ## Transformed additive genetic values - variant 1
## for(k in 1:nI) {
## u[k] ~ dnorm(0, tau2a)
## a[id[k]] <- u[k] * wsqr[id[k]] +
## 0.5 * (a[fid[k]] + a[mid[k]])
## }
## a[nU] <- 0 # NULL (zero) holder
##
## ## Location priors
## for(j in 1:nB) {
## b[j] ~ dnorm(0, 1.0E-6)
## }
##
## ## Phenotypes
## for(i in 1:nY) {
## y[i] ~ dnorm(mu[i], tau2e)
## mu[i] <- inprod(x[i, ], b[]) + a[idy[i]]
## }
## }
dataCholVar1 <- do.call('rbind',
lapply(c("./Results/Large_Example/ESSCholVar1.RDS",
"./Results/Large_Example/ESSCholVar1Var.RDS"), readRDS))
dataCholVar1$model <- as.factor(dataCholVar1$model)
g1 <- dataCholVar1 %>%
ggplot(aes(y= cost_mean, x= Generation)) +
geom_point(aes(colour = par)) +
geom_line(aes(linetype = par, colour = par)) +
xlab("Generation") +
ylab("Mean computing cost per sample (seconds)")+
ylim(0, max(dataChol$cost_max)) +
facet_wrap(~model)+
theme_bw() +
theme(legend.position = "none")
g2 <- dataCholVar1 %>%
ggplot(aes(y= cost_max, x= Generation)) +
geom_point(aes(colour = par)) +
geom_line(aes(linetype = par, colour = par)) +
xlab("Generation") +
ylim(0, max(dataChol$cost_max)) +
ylab("Maximum computing cost per sample (seconds)")+
facet_wrap(~model)+
theme_bw()
g1+g2
- Animal model variant 2 using Cholesky decompostion of the NRM
# Variance unknown
AnimalModelBUGSVariant2Var
## function()
## {
## ## Priors - precision parameters
## tau2e ~ dgamma(0.001, 0.001)
## tau2a ~ dgamma(0.001, 0.001)
## sigma2e <- pow(tau2e, -1)
## sigma2a <- pow(tau2a, -1)
##
## ## Transformed additive genetic values - variant 2
## sigmaa <- pow(tau2a, -2)
## for(k in 1:nI) {
## u[k] ~ dnorm(0, 1)
## a[id[k]] <- u[k] * wsqr[id[k]] * sigmaa +
## 0.5 * (a[fid[k]] + a[mid[k]])
## }
## a[nU] <- 0 # NULL (zero) holder
##
## ## Location priors
## for(j in 1:nB) {
## b[j] ~ dnorm(0, 1.0E-6)
## }
##
## ## Phenotypes
## for(i in 1:nY) {
## y[i] ~ dnorm(mu[i], tau2e)
## mu[i] <- inprod(x[i, ], b[]) + a[idy[i]]
## }
## }
dataCholVar2 <- do.call('rbind',
lapply(c("./Results/Large_Example/ESSCholVar2.RDS",
"./Results/Large_Example/ESSCholVar2Var.RDS"), readRDS))
dataCholVar2$model <- as.factor(dataCholVar2$model)
g1 <- dataCholVar2 %>%
ggplot(aes(y= cost_mean, x= Generation)) +
geom_point(aes(colour = par)) +
geom_line(aes(linetype = par, colour = par)) +
xlab("Generation") +
ylab("Mean computing cost per sample (seconds)") +
# ylim(0, max(dataCholVar2$cost_max)) +
facet_wrap(~model, scales = "free") +
theme_bw() +
theme(legend.position = "none")
g2 <- dataCholVar2 %>%
ggplot(aes(y= cost_max, x= Generation)) +
geom_point(aes(colour = par)) +
geom_line(aes(linetype = par, colour = par)) +
xlab("Generation") +
# ylim(0, max(dataCholVar2$cost_max)) +
ylab("Maximum computing cost per sample (seconds)") +
facet_wrap(~model, scales = "free") +
theme_bw()
g1+g2
Here we show that ‘prior’ Cholesky decomposition can solve a system of equations as an alternative to direct matrix inversion. However, it has an additional cost per effective sample with a maximum cost per sample of approximate \(82\) seconds against \(12.5\) seconds observed to direct matrix inversion. Cholesky decomposition is computationally costlier as the size of matrix grows due to its floating-point math operations per \([N \times N]\) matrix, which is generally estimated as \(\mbox{FLOPS}=4\frac{n^3}{3} = O(n^3)\) (Chari and Salon 2000). On the other hand, this approach can gain the inverse matrix more numeric efficiently than from the whole matrix as we only need to invert \(\boldsymbol{L}\), not \(\boldsymbol{A}=\boldsymbol{L}\boldsymbol{L}^{T}\) (Alfriend et al. 2010). Cholesky decomposition is also helpful for efficient numerical solutions. Its computational rate and efficiency highly depend on implementation details and the computing device’s architecture details used to run the analysis (Chari and Salon 2000).
We also observed that the animal model variant 2 using Cholesky decomposition of the NRM and Gibbs sampling algorithm has the highest mean and maximum cost per sample. Thus we recommend avoiding such parameterization to run MCMC with Gibbs sampling sampler. Our results showed the maximum cost per effective sample was 4400 seconds and 249 seconds for unknown and known variance components, respectively. Therefore, the original implementation of Damgaard ((Damgaard 2007b); variant 2 in function Variantes 2) is the most time-consuming approach, and the implementation with the transformed additive genetic values constantly produced chains with a considerably lower effective sample size when using Gibbs sampling sampler. Furthermore, Shariati and Sorensen (2009b) performed extensive tests of various McMC samplers for animal model and concluded that none of the samplers was the best in all scenarios, but transformed and blocked samplers could be beneficial when there are large posterior correlations between the parameters.